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Abstract 

A model of multicellular systems with several types of cells is developed from the phase field model. The 
model is presented as a set of partial differential equations of the field variables, each of which expresses 
the shape of one cell. The dynamics of each cell is based on the criteria for minimizing the surface area 
and retaining a certain volume. The effects of cell adhesion and excluded volume are also taken into 
account. The proposed model can be used to find the position of the membrane and/or the cortex of 
each cell without the need to adopt extra variables. This model is suitable for numerical simulations of a 
system having a large number of cells. The two-dimensional results of cell adhesion, rearrangement of a 
cell cluster, and chemotaxis as well as the three-dimensional results of cell clusters on the substrate are 
presented. 

Introduction 

In order to investigate the structural patterns of cellular systems, several cell models have been reported, 
including the vertex dynamics model [HIS], the center dynamics model OH], and the cellular Potts 
model [51IB| . Both the vertex dynamics model and the center dynamics model express cell patterns using 
polygons. In the vertex dynamics model, a cell or a cluster of cells is represented by a polygon formed 
by linking several vertices. Each vertex is driven by forces acting on it. This model has been adopted 
for morphogenesis in Xenopus notochords as well as cell deformation and rearrangement by applying 
mechanical forces [UlT]. In the center dynamics model, a node represents a cluster of cells and receives 
forces from its neighboring nodes. Cell aggregation, locomotion, rearrangement, and morphogenesis in 
vertebrate limb buds have been investigated using this model [3ll4ll8Hl0] . Although the mechanical pro- 
cesses during tissue developments can be well investigated, artificial treatments are required for numerical 
simulations in these models based on polygons. For example, in the vertex dynamics model, cell rear- 
rangement is realized by manually exchanging two vertices that approach each other jl]. In the center 
dynamics model, in order to express the cell division, it is necessary to add a new node in the vicinity of 
the existing node [4lll0j. 

In contrast, the cellular Potts model represents each cell as a cluster of grid points under the constraint 
of constant volume. Thus, the artificial treatments mentioned above arc not required for simulations in 
this model. We can investigate the deformation of an individual cell in a multicellular system using this 
model, considering the effects of excluded volumes and adhesions of the cells. This model successfully 
described several biological behaviors jllj . For example, numerical calculations with regard to cell sorting, 
biofilm formation, and chemotactic movement have been performed [5ll6l ll2[[T3] . However, running the 
simulations requires fluctuations, and the forces between cells are not expressed directly in this model. 

Therefore, we consider a new type of a model for multicellular systems, which is based on the phase 
field model. The effects of cell adhesion and excluded volume are taken into account. In the proposed 
model, the free energy is described in terms of a vector variable, the number of components of which is 
equivalent to the total number of cells in the system. The shape of one cell is expressed by one component 
of the vector variable. The time evolutions are described by a set of partial differential equations that are 
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obtained by taking the functional derivative of the free energy. Thus, fluctuations are not required for 
numerical simulations. In addition, by adopting auxiliary variables that arc used for calculation of the 
interactions between the cells, a program that consumes little computational memory can be designed. 
That is to say, the proposed model can be used to describe a system containing a large number of cells. 
The proposed model differs from previous models of multicellular systems in that the position of the cell 
membrane and/or cortex can also be expressed without the need to adopt extra variables because the 
phase boundary interface is treated as a diffuse interface of finite width using the phase field method. 

The phase field model has been applied to a wide range of problems, such as crystal growth [T4Ul8j . 
Very recently, the cell shape of the fish keratocyte has been modeled using this method, where the 
membrane bending force and the surface tension of the cell were considered [19]. However, to our 
knowledge, this is the first report applying the phase field method to the multicellular system. 



Results 
Model Equation 

We consider a multicellular system containing several types of cells and allow changes in the size and 
adhesive strength of each cell type. As a first step, we express the shape of one cell using the phase field 
method. 

The following Ginzburg-Landau free energy is considered: 



E[u] = I 



dv+^{Vo-v{u))\ (1) 



where 17 denotes the area of the system, and the coefficients Dqj Q^Oj and Vb are positive constants. The 
variable u{r,t) is an order parameter referred to as the phase field, where r is the position, and t is the 
time. The function v{u) is given as 



v{u) ^ / h{u)dr, (2) 



where the function h{u) is defined as 



h{u) =u\3-2u). (3) 



By taking the functional derivative of Equation [T] with respect to u, the time evolution of u is derived as 
follows: 

du _ SE 
dt Su 

= Do\/^u + u{l-u)(u~^ + fo{u)y (4) 
fo{u) = ao{Vo-v{u)), (5) 

where t is a positive constant. Equation |4] guarantees the monotonic decrease in the free energy. 

If the function /o is set to be constant /, then Equation |4] is referred to as the Allen-Cahn equation 
in the field of materials science and is known for having a smooth front solution connecting the regions 
u = I and u = 0. The AUen-Cahn equation can easily be solved in one dimension as u = {1 — tanh[(a; — 
Vt) / {2y/2Do)]} /2, where the front velocity V = \f2J\j jr. This means that the front moves such that 
the region of u = 1 (u = 0) expands if / > (/ < 0). 

Note that the function v{u) can be regarded as the volume of the region in which u = 1 because 
h{\) = 1 and /i(0) = 0. Therefore, as discussed above. Equation [S] indicates that the region of m = 1 
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expands (shrinks) until v{u) ~ Vo when v(u) < Vq {v(u) > Vo). The last term of Equation [U which has 
a minimum at v{u) = Vq, also expresses the constraint of the constant volume of u. 

As shown in Figure [TJ the region of m = 1 takes the form of a circle in two dimensions and a sphere in 
three dimensions in the steady state. Thus, the shape of the cell in the simplest case can be described by 
a single-order parameter u, such that u > Uceii in the region with the cell (< Uceii in the region not taken 
up by the cell) with a constant u^eu G (Oj !)■ Based on the fact that u has an interface with a thickness 
on the order of \/Dq, the cell cortex can also be expressed as a function of u, e.g.. u(l — u) (see Figure 

In order to describe the multicellular system, a vector variable u(r,t) = (ui(r,i),--- ,UA/(r,i)) is 
considered, where M is the total number of cells in the system. The component Mm(r, t) (to = 1, • • • , M) 
describes the shape of the TO-th cell. We also use the variable s(r,t) to represent the shape of substances 
interacting with the cells, such as the wall (Figures [5] and [7]), the substrate (Figures H]) , and the ECM. 

The model free energy for the multicellular system is written as 

^[u, s] = E,M [u] + E,nt [u] + Es [u, s] , (6) 

where Eceii determines the shape of the cell, Eint describes the interactions between each cell, and Eg 
expresses the interactions between the cells and substances external to them. The form of E^eii is obtained 
by modifying Equation [T] using the vector variable as follows: 



Ecell [u] = 



E 



E 



D{im) |2 , 1 _ 2 



dr 



12 



-{V{l,n)-V{u^))\ (7) 



where £m is the cell type of the w-th cell. The coefhcients D{(), a(£), and V{t} {i = 1, ■ ■ ■ L) are positive 
constants, where L is the total number of cell types in the system. As discussed in the paragraph below 
Equation [51 Equation [7] indicates that the thickness of the cell interface is on the order of ^jD{li) and 
that the speed at which the volumes of the type-£ cells approach the target volume V{t} is controlled by 
the value of a{t}. That means a{t) determines the cell size growth. Here, E^nt can be presented in the 
following form: 



Eint[v] = E E ^^^"'q^™ h{ujn)h{um')dr 



EE 



E 



7(0 
12 



|V/i(u™)|2dr, (8) 



where (3{£,£'), '/]{£,£'), and 7(f) {£,£' = 1, - ■ ■ L) are positive constants. The first term on the right-hand 
side of Equation [5] represents the effect of the excluded volume by increasing the energy if the cells 
overlap, whereas the second term represents the effect of cell adhesion by decreasing the energy if the 
cell cortices overlap. This adhesion term becomes negative in the region in which cell adhesion occurs. 
In order to prevent divergence due to this adhesion term, we introduce the third term on the right-hand 
side of Equation [8] with the condition whereby j{£) > ri{£, £). Similarly, the interaction between cells and 
substances external to the cells is expressed as follows: 



Es[n,s\ = / h{u^)h{s)dv 



+ E / ^^("™) • '^K^)dr. (9) 
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where /3s(€) and ris{£) = 1, • • • ,L) are positive constants. 

Taking the functional derivative of Equation [5] with respect to u„i, the fohowing time evolution 
equations are obtained: 



T,, 



dt 



2„ 



+ 9int{u 

), (10) 

f{ujn,s,(f>) = a{e„i)(y{£„i) - v{um)) 

- ^/3(^™,£)[0£-/lK„)5£„,£]-/3s(OM5), (11) 
gint{Um,(t>) = ^?7(^m,-^)V - U„)V{0£ - /l(u„)(5f„,f}] 

+ 7(C)V[u™(l-M„,)V/i(u„)], (12) 

5s(Um,s) = r]s{£m)y [Um{l - Urn)yh{s)] , (13) 

where t„ is a positive constant, and is the Kronecker delta, which is 6ij = 1 {Si^j = 0) if i = j (i 7^ j). 
The vector variable 4'{r,t) = ((/)i(r,<), • • • ,(/)L(r,<)) is an auxiliary variable that is defined as follows: 

t)=J2 Kum{v, t))6i^,i. (14) 

ni 

As shown in Figure [21 the region occupied by the type--^ cells can be identified by <f>i. 

Note that the interaction terms in Equation [10] are not written explicitly in terms of the variables 
Um' {rn' = 1, • • • , M ^ m) but are instead written in terms of the auxiliary variable (j). Moreover, the 
components Eint and can also be presented in terms of </>, as follows: 



— / h{um) dv 

rn •'^ 

PI 'J ^1 



/ |V/i(u™)|^dr, (15) 



+ Y.Hr f ^'^(')dr. (16) 
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We adopted the second term on the right-hand side of Equation [TT] and the first term on the right-hand 
side of Equation [12] to express the excluded volumes and the cell adhesions, respectively, because these 
terms are the simplest among the several alternatives, which can be written in terms of (j) both in the 
time evolution equation for u and the component Eint- 
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Numerical implementation 

In order to rapidly simulate a system having numerous cells, it is important to design a program that 
does not consume a large amount of computational memory and to increase the simulation speed. These 
two requirements are easily satisfied because Equation [10] is not written explicitly in terms of Um' ("z' = 
1, • • • , A/ 7^ m). Once (j) is obtained for each time step, the time evolution of can be computed 
independent of Um'- Such a program is fully compatible with parallel computation. Moreover, the shape 
of the TO-th cell can be obtained by computing the equation for within the small region il,„, which 
covers the region of w„i > 0. This reduces the computational memory and increases the simulation speed. 
The position r,„(t), which indicates the center position of 17„i measured for the entire system, must be 
moved along with the movement of the center position of the m-th cell. Since Um = is realized outside 
the region 17™, the Dirichlet boundary condition must always be set for the small region fim. 

The estimation of the required memory is described below. Since the number of cell types L is generally 
much smaller than the number of cells Af, the memory increase by introducing </> becomes smaller than 
the memory decrease by computing m,„ within the small region ^Im- For simplicity, we assume that each 
of the cells has the same volume, i.e., V{1) = ■ ■ ■ = V{L) = V and that the entire system is covered by 
the cells, i.e., ~ VM. Then, the computational memories for u, r{t) = (ri(i), • • • ,rj^j{t)), and (f> are 
roughly estimated as dM, and LVM/5'^, respectively, where d is the spatial dimension and 5 

is the size of the spatial grid. Therefore, the total memory required to compute Equation 1101 using (j) is 
linearly dependent on M. On the other hand, in order to compute cell-cell interactions without using 
(j), the value of u must be preserved over the entire region f2. Then, the computational memories for 
solving Equation [TOl increase by VM/6''' x M oc M^. These estimations reveal that the introduction of (j) 
is very useful for computation in the case of a system that contains a large number of cells, even in three 
dimensions. 

Numerical Simulation 

Figure 13] shows the numerical results for two cells of the same type, i.e., M = 2 and ii = £2 ^ L — 1, 
with different adhesion strengths. The curves in the top row of the graphs indicate the contour lines 
of Ujn = 0.2 (m = 1 and 2), and the x symbols indicate the positions of the centers of the cells = 
(/j^ YUmdr) I (Jj^j Urndr). The variable e,,(r, t) is given as follows: 

e„(r, i) = E E ^^^^^^ V;z(^„,(r, t)) ■ ^h[u„, (r, t)). (17) 

The integral over r of is identical to the second term on the right-hand side of Equation [S] The u,„ 
and e,, profiles along the dotted line in the top row have been plotted in the middle and bottom rows of 
the graphs, respectively. Since has a non-zero value only in regions in which cell adherence occurs, e,, 
is an indicator of locations at which cell adherence occurs. Initially, the distance between the centers of 
cells is set to 1.6000. After a sufficiently long simulation time (i = 50,000), the two cells move closer to 
each other as the value of 77(1, 1) increases, such that the distances between the cell centers are 1.662 in 
the case of Panel A with 77(1, 1) = 0.0000, 1.355 in the case of Panel B with 77(1, 1) = 0.004, and 1.156 in 
the case of Panel C with 77(1, 1) = 0.008. 

Figure [4] shows snapshots of three-dimensional simulations at i = 500. The periodic boundary 
conditions are imposed on Vt. The solid substrate is introduced by setting the variable s as s(r) = 
(1 — tanh((2; — Zf) / ef)) /2, where Zf and e/ are positive constants. The light gray surfaces are contour 
plots of Um — 0.1 [m — 1, - ■ ■ , 10) and the dark gray surfaces represent contour plots of s = 0.1. We set 
77(1, 1) as 0.0000, 0.0100, and 0.0219 for the simulations shown in Panels A, B, and C, respectively, where 
the other parameters are the same for all cases. If the cell adhesions are weak, the cells push against each 
other, and their positions are determined as shown in Panel A. On the other hand, for the case in which 
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the cell adhesions are sufficiently strong, the cell positions are decided by the pulling force between cells, 
and the surface of the cell layer becomes flat, as shown in Panel C. 

Figure [5] shows the numerical results for cell deformation and rearrangement. A cell cluster oi AI — 8 
and L = 2 is sandwiched between two walls that move at a constant speed. In this calculation, considering 
the variable s as an order parameter that corresponds to the walls, the time evolution of s is calculated 
as s{r,t) = 1 — (1 + tanh((x — xi{t))/es)) (1 — tanh((a; — Xr{t))/es)) /4, where Eg is a positive constant. 
The locations of the left and right walls are denoted as xi{t) and Xr{t), respectively. Panel A shows the 
results for the case in which the adhesion strength between cells of the same type is stronger than that 
between cells of different types {ri{l, 1) = ri{2, 2) ~ 0.008 and 2) = 0.005), whereas Panel B shows the 
results for the opposite ease (77(1, 1) = 77(2,2) = 0.005 and 77(1,2) = 0.008). Light gray, dark gray, and 
black areas represent the positions of the type-1 cells, the type-2 cells, and the walls, respectively. Cells 
adhering to the walls are stretched by the moving walls, causing the cells to be deformed and rearranged. 
Cells that are rearranged as weakly adhered cells detach first. In Panel B, the cells separate into two 
groups at approximately t = 22,000 and relax to almost their original shape at t = 26,000. The time 
evolution of the total energy E is plotted in Figure [6] The solid line shows the results for Panel A of 
Figure [5l and the dotted line shows the results for Panel B of Figure [5] There is no monotonic decrease 
in total energy because the walls stretch the cell clusters. Comparison of Figures [5] and |6] reveals that the 
energy decreases significantly when cell rearrangement occurs. 

Finally, we show that the additional cell behavior can also be incorporated into the proposed model. 
For example, the chemotactic movement of the cell can be described by adding a new term, such as 
gchem = — M(^m) V • («,„ Vc) to the right-hand side of EguationfTOl where the variable c(r, t) is the chemical 
concentration in extracellular regions. The parameter ii{£m) indicates the sensitivity of the m-th cell to the 
gradient of c. Figure [7] shows the time evolution of a system with cells having chemotaxis. Light gray and 
dark gray represent type-1 and type-2 cells, respectively. In this case, we consider the variable s as an order 
parameter that corresponds to the wall. The fifty cells are surrounded by the unmoving wall defined as 
s(r) = 2-(l-Ftanh((x-a;i)/eu,))(l-tanh((a;-Xr)/e^))/4-(l-Htanh((y-yf,)/e„))(l-tanh((?/-yt)/e,„))/4, 
where Xi, Xr, yt, Ut, and are positive constants. Cell adhesion is not considered in this simulation. By 
setting /^(l) = 0.0 and /i(2) — 1.0, it is assumed that the only type-2 cells can sense the gradient of the 
chemical concentration c. For simplicity, the form of c is assumed not to be affected by Um or t and is 
taken as c = cqx, where cp is a constant. It is found numerically that type-2 cells move toward the c-rich 
region by pressing against type-1 cells. 

Discussion 

We proposed a new type of cell model based on a phase field model, including the effects of excluded 
volumes and cell adhesions. The proposed model is based on a concept similar to the cellular Potts 
model, but the time evolutions of cell shapes in the proposed model differ from those in the cellular 
Potts model. In the cellular Potts model, the time evolutions of the spins are computed by the Monte 
Carlo method, and thus the fluctuations are required for the time evolution. On the other hand, the 
time evolution equations in the present model are written in the form of partial differential equations, 
whereby fluctuations are not necessary in order to run the simulations. In addition, the proposed model 
is thought to be more appropriate for investigating problems in which a small volume variant must be 
accounted for, because the proposed model is continuous in any parameter. 

Since the cell shapes are represented by interfaces of finite thickness, the proposed model has the 
potential to be applied to the investigation of not only shape changes due to interactions between cells 
(Figures [3] and HI) and rearrangements of cells in clusters (Figure [S|) but also phenomena requiring knowl- 
edge of the position of the cell membrane and/or cortex. It is easy to incorporate additional cell behaviors 
such as chemotaxis (Figure [7]) into the proposed model by adding corresponding terms. At the stage of 
numerical implementation, a program that is suitable for parallel computing and that consumes little 
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computational memory can be designed by introducing the auxiliary variable 0, which is commonly used 
for the calculation of interactions between cells. Therefore, simulations can be performed for a system 
with numerous cells, even in three dimensions (Figure U). 

The proposed model can express the time evolution of changes in cell shape due to the interactions 
between cells, cell differentiation by changing the cell type, cell size growth, cell movement, and cell 
death by deleting the corresponding component of u. Thus, this model may well provide a useful tool 
for approaching the problem of morphogenesis, although this remains a subject for future study in order 
to estimate the parameters by comparison with earlier models and with experimental data. We plan to 
include the cell division, in the process of which the cortex of the dividing cell is known to be important 
[20H22] . in the present model and to approach the problem of morphogenesis. 
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Figure Legends 
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Figure 1. Shape of the phase field u. The integral of u over r is set to be maintained. Panel A; top 
view. Panels B and C: profiles of u and u(l — u) at the centerline in Panel A, respectively. 
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Figure 2. Schematic diagram of u and (p. In Panel A, type-1 {m ~ 1 and 2) and type-2 (to = 3 
and 4) cells are represented by gray and black circles, respectively. The contours of 0i and 02 are 
indicated by curved lines in Panels B and C, respectively. 
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Figure 3. Two-dimensional results of cell adhesions. The case of two cells {M = 2) of the same 
type {L ~ 1) is considered. Numerical calculations were performed with 77(1, 1) = 0.000 in Panel A, 
77(1, 1) = 0.004 in Panel B, and 77(1, 1) = 0.008 in Panel C. The top row shows contour plots of Um — 0.2 
(to = 1,2). The X symbol indicates the centers of cells. The middle and bottom rows show the profiles 
of Um and e,, along the dotted line shown in the top row. The size of the simulation box is f2 = 5 x 5, 
and the size of the spatial grid is (5 0.05. The time increment is dt = 0.01. The remaining parameters 
are set as follows: t„ = 1, £1(1) = 0.001, V{1) = 1, a(l) = 1, 7(1) = 0.01, ;3(1, 1) = 1, and 
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Figure 4. Three-dimensional results of cell adhesions on the substrate. The case of 10 cells 
{M ~ 10) of the same type (L = 1) is considered. Numerical calculations were performed with 
77(1, 1) 0.0000 in Panel A, 77(1, 1) = 0.0100 in Panel B, and 77(1, 1) = 0.0219 in Panel C. Light and 
dark gray surfaces are contour plots of Um = 0.1 {m = 1, • • • , 10) and s = 0.1, respectively. The 
diagonal, top, and side views for each result are shown in the top, middle, and bottom rows, 
respectively. The size of the simulation box is 17 = 5 x 5 x 4, and the size of the spatial grid is (5 = 0.05. 
The time increment is dt = 0.01. The remaining parameters are set as follows: r„ = 1, D{1) ~ 0.001, 
V{1) = 2.26, a(l) = 100, 7(1) = 0.022, /J(l, 1) = Psil) = 1, ?7s(l) = 0.01, Zf = 0.8, and e/ = 2Vomi. 



12 





t=7600 



t=6000 




t=20000 



t=20000 




t=26000 



t=20800 




t=34000 



t=26000 



Figure 5. Two-dimensional results of cell deformation and rearrangement in a cluster. The 

cluster is composed of eight cells (M — 8) of two types (L = 2). Light and dark gray areas represent the 
region of Um > 0.2. Light gray areas indicate the locations of type-1 cells, and dark gray areas indicate 
the locations of type-2 cells. Black areas represent the walls (s > 0.5). Numerical calculations were 
performed with r]{l, 1) = 77(2, 2) = 0.008 and 77(1, 2) = 0.005 in Panel A and 77(1, 1) = r/(2, 2) = 0.005 
and 77(1, 2) = 0.008 in Panel B. The left and right walls are assumed to move at a uniform velocity, 
xi — 7 — Vgt, Xr = 13 + Vgt, and = 0.0001. The size of the simulation box is = 20 x 15, and the 
size of the spatial grid is (5 = 0.05. The time increment is dt = 0.01. The remaining parameters are set 
as follows: Tu = 1, Du{l) = Du{2) = 0.001, y(l) = V{2) = 4, a(l) = a{2) = 10, 7( 1) =7 (2) = 0.01, 
/3(1, 1) = /3(1, 2) = /3(2, 2) = 13,(1) = /3,{2) = 0.1, 77,(1) = 77,(2) - 0.01, and = 2^0:002 . 



0.40 



0.39 



0.38 



0.37 



0.36 







A 

B 




m/ ^A, 











10000 20000 30000 

i 



Figure 6. Plots of the total energy E with respect to time. The solid line shows the results for 
Figure [IK, and the dotted line shows the results for Figure [33. 
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Figure 7. Two-dimensional results of chemotactic movement of cells. The case of fifty cells 
(A/ = 50) of two types (L = 2) is considered. Light gray (dark gray) areas indicate the region of 
M„i > 0.2, for the case in which the m-th cell is a type-1 (type-2) cell. Black areas represent the walls 
(s > 0.5). Numerical calculation was performed with ^(1) = 0.0 and /i(2) = 1.0. The other parameters 
are set as follows: size of the simulation box = 10 x 10, size of the spatial grid 8 = 0.05, time 
increment dt = 0.01, t„ = 1, £)„(1) = D^{2) = 0.001, 1/(1) = V{2) = 1, a(l) = a(2) = 1, 
/?(1, 1) = /3(1, 2) - /3(2, 2) = /3,(1) = /3,(2) = 1, 7(1) = 7(2) = 0, 

?7(1, 1) = r/(l, 2) = 7/(2, 2) = r/,(l) = r/,,(2) = 0, = 270002, cq - 0.01, xi=y^,^ 0.8, and 
Xr = yt ^ 9.2. 



